
// Table 1: Effects of road tolls and bus lane distance on household BEV ownership, 2015-2017

*=========================================
// Step 1: Run regressions and save estimates
*=========================================

// Loop over couple (1) and single-adult (0) HH
forvalues i=0(1)1 {

use    "${dataout}MainDataset" , clear
keep if year>=2015 & year<=2017
keep if couple==`i'
keep  $KeepOutcome  $KeepTreat  $KeepControl 

	if `i'==0 {
		//* Single-adult: Generate one value per row values 
		egen employed = rowmax(employed1 employed2) 
		egen retired = rowmax(retired1 retired2) 
		egen age = rowmax(age1 age2) 
		egen wies_decile = rowmax(wies1_decile wies2_decile) 
		egen wealth_decile = rowmax(wealth1_decile wealth2_decile) 
		egen grputd = rowmax(grputd1 grputd2) 
		egen time = rowmax(time1 time2) 
		egen dist = rowmax(dist1 dist2) 
		egen grk_bed_num  = rowmax(grk_bed1_num  grk_bed2_num) 
		egen PublicTransitTime  = rowmax(PublicTransitTime_1  PublicTransitTime_2) 
		egen VENTETID  = rowmax(VENTETID_1  VENTETID_2) 
		egen BOARDINGS  = rowmax(BOARDINGS_1  BOARDINGS_2) 
		egen TAKST  = rowmax(TAKST_1  TAKST_2 )
		
		local household     $householdS
		local employment    $employmentS 
		local age 		    $ageS
		local income 	    $incomeS 
		local education	    $educationS
		local distance 	    $distanceS   
		local time 		    $timeS  
		local FE		    $FES 
		local publictime	$publictimeS   
		local publicquality $publicqualityS
		local clustervar	$clustervarS
	}
	
	if `i'==1 {
		
		local household     $household
		local employment    $employment
		local age 		    $age
		local income 	    $income 
		local education	    $education
		local distance 	    $distance   
		local time 		    $time 
		local FE		    $FE
		local publictime	$publictime   
		local publicquality $publicquality
		local clustervar	$clustervar
	}


** Select variables to loop over: 
local outcomelist BEV_fam_yes  
local treatvar1list toll_fam_mean
local treatvar2list ptl_fam_km_mean     


** Regressions - loop over: outcome, treatment, year
foreach outcome in `outcomelist'{
foreach treatvar1 in `treatvar1list'{
foreach treatvar2 in `treatvar2list'{

// g4
reghdfe `outcome' `treatvar1'  `treatvar2'  `age'    `distance'  `time'  `publictime'  `publicquality'   , ///
absorb(`FE'   `household'  `income'  `employment'  `education'  i.year   ) vce(cluster `clustervar') 
*
summarize `outcome' if e(sample)==1
estadd scalar MeanDep = r(mean), replace
summarize `treatvar1' if e(sample)==1
estadd scalar Treat1 = r(mean), replace
summarize `treatvar2' if e(sample)==1
estadd scalar Treat2 = r(mean), replace
*
estadd local year="2015-17" , replace
estadd local grkFE="\checkmark", replace
estadd local grkbFE="\checkmark", replace
estadd local HHcon="\checkmark", replace
estadd local HHWork="\checkmark", replace
estadd local PublicTime="\checkmark", replace
estadd local grkFE_W="Yes", replace
estadd local grkbFE_W="Yes", replace
estadd local HHcon_W="Yes", replace
estadd local HHWork_W="Yes", replace
estadd local PublicTime_W="Yes", replace
*
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estimates save "${ster}reg_g4_2015_2017_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple`i'" , replace


** Mark main sample
capt drop MainSample
gen MainSample=0
replace MainSample=1 if e(sample)==1

// g1
reghdfe `outcome' `treatvar1'  `treatvar2'   if MainSample==1 , absorb(`FE'   i.year  ) vce(cluster `clustervar') 
*
summarize `outcome' if e(sample)==1
estadd scalar MeanDep = r(mean), replace
summarize `treatvar1' if e(sample)==1
estadd scalar Treat1 = r(mean), replace
summarize `treatvar2' if e(sample)==1
estadd scalar Treat2 = r(mean), replace
*
estadd local year="2015-17" , replace
estadd local grkFE="\checkmark", replace
estadd local grkbFE="\checkmark", replace
estadd local HHcon="", replace
estadd local HHWork="", replace
estadd local PublicTime="", replace
estadd local grkFE_W="Yes", replace
estadd local grkbFE_W="Yes", replace
estadd local HHcon_W="", replace
estadd local HHWork_W="", replace
estadd local PublicTime_W="", replace
*
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
*
estimates save "${ster}reg_g1_2015_2017_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple`i'" , replace


// g2
reghdfe `outcome' `treatvar1'  `treatvar2'  `distance'  `time'   if MainSample==1 , ///
absorb(`FE'  i.year  ) vce(cluster `clustervar') 
*
summarize `outcome' if e(sample)==1
estadd scalar MeanDep = r(mean), replace
summarize `treatvar1' if e(sample)==1
estadd scalar Treat1 = r(mean), replace
summarize `treatvar2' if e(sample)==1
estadd scalar Treat2 = r(mean), replace
*
estadd local year="2015-17" , replace
estadd local grkFE="\checkmark", replace
estadd local grkbFE="\checkmark", replace
estadd local HHcon="", replace
estadd local HHWork="\checkmark", replace
estadd local PublicTime="", replace
estadd local grkFE_W="Yes", replace
estadd local grkbFE_W="Yes", replace
estadd local HHcon_W="", replace
estadd local HHWork_W="Yes", replace
estadd local PublicTime_W="", replace
*
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estimates save "${ster}reg_g2_2015_2017_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple`i'" , replace

}
}
}
}


*=========================================
// Step 2: Load estimates saved to disk + make calculations
*=========================================

* Load data 

use    "${dataout}MainDataset" , clear
keep if year >= 2015

*----- Used for single-adult hh ---------
egen employed = rowmax(employed1 employed2) 
egen retired = rowmax(retired1 retired2) 
egen age = rowmax(age1 age2) 
egen wies_decile = rowmax(wies1_decile wies2_decile) 
egen wealth_decile = rowmax(wealth1_decile wealth2_decile) 
egen grputd = rowmax(grputd1 grputd2) 
egen time = rowmax(time1 time2) 
egen dist = rowmax(dist1 dist2) 
egen grk_bed_num  = rowmax(grk_bed1_num  grk_bed2_num) 
egen PublicTransitTime  = rowmax(PublicTransitTime_1  PublicTransitTime_2) 
egen VENTETID  = rowmax(VENTETID_1  VENTETID_2) 
egen BOARDINGS  = rowmax(BOARDINGS_1  BOARDINGS_2) 
egen TAKST  = rowmax(TAKST_1  TAKST_2) 
*---------------------------------------------------------

local yearlist  2015_2017   
local outcomelist   BEV_fam_yes 
local treatvar1list toll_fam_mean     
local treatvar2list ptl_fam_km_mean     


capt eststo drop *
capt estimates drop *

foreach outcome in `outcomelist'{
foreach treatvar1 in `treatvar1list'{
foreach treatvar2 in `treatvar2list'{
foreach y in `yearlist' {

// g1 - couples
eststo drop *

estimates use "${ster}reg_g1_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple1" 
* Load in approximately correct sample
estimates esample: `outcome' `treatvar1' `treatvar2' year children antpers_i_regstat_famnr secondhome employed retired age wies_decile  wealth_decile grputd dist time grkrets_num grk_bed_num PublicTransitTime VENTETID BOARDINGS  TAKST if couple == 1, replace 
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
eststo

// g2 - couples
estimates use "${ster}reg_g2_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple1" 
* Load in approximately correct sample
estimates esample: `outcome' `treatvar1' `treatvar2' year children antpers_i_regstat_famnr secondhome employed retired age wies_decile  wealth_decile grputd dist time grkrets_num grk_bed_num PublicTransitTime VENTETID BOARDINGS  TAKST if couple == 1, replace 
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd local HHcon_W="", replace   /*Fix a mistake..*/
eststo

// g4 - couples
estimates use "${ster}reg_g4_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple1" 
* Load in approximately correct sample
estimates esample: `outcome' `treatvar1' `treatvar2' year children antpers_i_regstat_famnr secondhome employed retired age wies_decile  wealth_decile grputd dist time grkrets_num grk_bed_num PublicTransitTime VENTETID BOARDINGS  TAKST if couple == 1, replace 
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
eststo

// g1 - single-adult
estimates use "${ster}reg_g1_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple0" 
estadd local year="2015-17" , replace
* Load in approximately correct sample
estimates esample: `outcome' `treatvar1' `treatvar2' year children antpers_i_regstat_famnr secondhome employed retired age wies_decile  wealth_decile grputd dist time grkrets_num grk_bed_num PublicTransitTime VENTETID BOARDINGS  TAKST if couple == 0, replace 
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
eststo

// g2 - single-adult
estimates use "${ster}reg_g2_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple0" 
* Load in approximately correct sample
estimates esample: `outcome' `treatvar1' `treatvar2' year children antpers_i_regstat_famnr secondhome employed retired age wies_decile  wealth_decile grputd dist time grkrets_num grk_bed_num PublicTransitTime VENTETID BOARDINGS  TAKST if couple == 0, replace 
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2)) 
estadd local Sample="Singles" , replace
estadd local HHcon_W="", replace   /*Fix a mistake..*/
estadd local year="2015-17" , replace
eststo

// g4 - single-adult
estimates use "${ster}reg_g4_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple0" 
estadd local year="2015-17" , replace
* Load in approximately correct sample
estimates esample: `outcome' `treatvar1' `treatvar2' year children antpers_i_regstat_famnr secondhome employed retired age wies_decile  wealth_decile grputd dist time grkrets_num grk_bed_num PublicTransitTime VENTETID BOARDINGS  TAKST if couple == 0, replace 
sum `outcome' if `treatvar1' == 0 & `treatvar2' == 0 & e(sample) == 1
estadd scalar rel_effect1 = _b[`treatvar1'] / r(mean)
estadd scalar rel_effect1se = _se[`treatvar1'] / r(mean)
estadd scalar rel_effect2 = _b[`treatvar2'] / r(mean)
estadd scalar rel_effect2se = _se[`treatvar2'] / r(mean)
estadd scalar elast1 = (_b[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast1se = (_se[`treatvar1'] /e(MeanDep)) / ( 1/e(Treat1))
estadd scalar elast2 = (_b[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2))
estadd scalar elast2se = (_se[`treatvar2'] /e(MeanDep)) / ( 1/e(Treat2)) 
eststo
}
}
}
}

est dir

*=========================================
// Step 3: Make Table 
*=========================================

local treatvar1 toll_fam_mean     
local treatvar2 ptl_fam_km_mean   

// Latex 
esttab est*  using  "${tables}Table1.tex" ///
, keep(`treatvar1'*  `treatvar2'*)   ///
title("Effects of road tolls and bus lane distance on household BEV ownership, 2015-2017.") ///
nonotes  nomtitle   /*label*/  coeflabels(toll_fam_mean "Road toll (NOK)" toll_fam_mean_KPI "Road toll (NOK)" toll_fam_sum_KPI  "Road toll, sum (NOK)"  ptl_fam_km_mean "Bus lane (km)") ///
booktabs   mgroups("Couple households" "Single-adult households" , pattern(1 0 0 1 0 0)   ///
prefix(\multicolumn{@span}{c}{)   suffix(}) span erepeat (\cmidrule(lr){@span})) ///
stats(N year grkFE grkbFE   HHWork HHcon PublicTime MeanDep Treat1 Treat2 midrule rel_effect1  elast1  rel_effect2  elast2 , label("N" "Year"  "Neighborhood residence FE" "Neighborhood work FE"  "Ind-work controls" "HH and ind controls"  "Public transit (minutes)"  "Mean outcome" "Mean road toll (NOK)" "Mean bus lane (km)" "\midrule \\ \noalign{\vspace{-10mm}}" "Road toll, prop. effect"  "Road toll elasticity"  "Bus lane, prop. effect"  "Bus lane elasticity")   fmt(%9.0fc a3)  ) ///
b(5) se(5) star(* 0.10 ** 0.05 *** 0.01) replace

// Latex - without table opening
esttab est*  using  "${tables}Table1_Fragment.tex" ///
, keep(`treatvar1'*  `treatvar2'*)   ///
nonotes  nomtitle   coeflabels(toll_fam_mean "Road toll (NOK)" toll_fam_mean_KPI "Road toll (NOK)" toll_fam_sum_KPI  "Road toll, sum (NOK)" ptl_fam_km_mean "Bus lane (km)") ///
booktabs   mgroups("Couple households" "Single-adult households" , pattern(1 0 0 1 0 0)   ///
prefix(\multicolumn{@span}{c}{)   suffix(}) span erepeat (\cmidrule(lr){@span})) ///
stats(N year grkFE grkbFE   HHWork HHcon PublicTime MeanDep Treat1 Treat2 midrule rel_effect1  elast1  rel_effect2  elast2 , label("N" "Year"  "Neighborhood residence FE" "Neighborhood work FE"  "Ind-work controls" "HH and ind controls"  "Public transit (minutes)"  "Mean outcome" "Mean road toll (NOK)" "Mean bus lane (km)" "\midrule \\ \noalign{\vspace{-10mm}}" "Road toll, prop. effect"  "Road toll elasticity"  "Bus lane, prop. effect"  "Bus lane elasticity")   fmt(%9.0fc a3)  ) ///
b(5) se(5) star(* 0.10 ** 0.05 *** 0.01) replace





